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Abstract. Analysis of the gravity and topography of Mars presently provides our primary 
quantitative constraints on the internal structure of Mars. We present an inversion of the long- 
wavelength (harmonic degree < 10) gravity and topography of Mars for lateral variations of 
mantle temperature and crustal thickness. Our formulation incorporates both viscous mantle 
flow (which most prior studies have neglected) and isostatically compensated density anomalies 
in the crust and lithosphere. Our nominal model has a 150-km-thick high-viscosity surface 
layer over an isoviscous mantle, with a core radius of 1840 km. It predicts lateral temperature 
variations of up to a few hundred degrees Kelvin relative to the mean mantle temperature, with 
high temperature under Tharsis and to a lesser extent under Elysium and cool temperatures 
elsewhere. Surprisingly, the model predicts crustal thinning beneath Tharsis. If correct, this 
implies that thinning of the crust by mantle shear stresses dominates over thickening of the 
crust by volcanism. The major impact basins (Hellas, Argyre, Isidis, Chryse, and Utopia) are 
regions of crustal thinning, as expected. Utopia is also predicted to be a region of hot mantle, 
which is hard to reconcile with the surface geology. An alternative model for Utopia treats it 
as a mascon basin. The Utopia gravity anomaly is consistent with the presence of a 1.2 to 1.6 
km thick layer of uncompensated basalt, in good agreement with geologic arguments about the 
amount of volcanic fill in this area. The mantle thermal structure is the dominant contributor to 
the observed geoid in our inversion. The mantle also dominates the topography at the longest 
wavelengths, but shorter wavelengths (harmonic degrees >4) are dominated by the crustal 
structure. Because of the uncertainty about the appropriate numerical values for some of the 
model’s input parameters, we have examined the sensitivity of the model results to the 
planetary structural model (core radius and core and mantle densities), the mantle’s viscosity 
stratification, and the mean crustal thickness. The model results are insensitive to the specific 
thickness or viscosity contrast of the high-viscosity surface layer and to the mean crustal 
thickness in the range 25 to 100 km. Models with a large core radius or with an upper mantle 
low-viscosity zone require implausibly large lateral variations in mantle temperature. 


Introduction 

In the absence of a lander network providing seismic and 
heat flow data, gravity and topography are our primary quanti- 
tative constraints on the internal structure of Mars. In this 
work, we present an inversion of Mars gravity and topography 
for mantle temperature anomalies and crustal thickness varia- 
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tions. We consider spherical harmonics up to degree 10, 
corresponding to wavelengths longer than 2100 km. Accord- 
ingly, our analysis emphasizes the major landforms on Mars: 
Tharsis, Elysium, Valles Marineris, and the large impact basins. 

Our model incorporates the effects of viscous mantle 
flow, something that has not been included in previous inver- 
sions of the Martian gravity field. Our approach is generally 
similar to the Venus model of Herrick and Phillips [1992], 
although our model differs from theirs in some details. We be- 
gin by describing the spherical harmonic data sets used in our 
inversion, along with a brief consideration of some of the sta- 
tistical properties of the harmonic models. We then describe 
the details of our model’s formalism. Next, we present results 
for a series of models, allowing for a plausible range of values 
for the core radius, mantle viscosity model, and mean crustal 
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thickness. We conclude with a brief consideration of the rela- 
tionship between the stress field predicted by our model and the 
observed tectonics of the Tharsis region. 


Data Sets 


The gravity field of Mars has been estimated using 
Doppler tracking of the Mariner 9 and Viking 1 and 2 space- 
crafts. In this work, we use the spherical harmonic degree 50 
model, GMM1, of Smith et al [1993]. The gravitational po- 
tential, U, for this field is 


f/(0,<(>,r) = 


GM 


50 / 

i + L I 

1 = 2 «i = 0 


R 

r 


P/,„(cos0) 


x ( C( m cos m<|> + S},„ sin m ( |> ) ( 1 ) 


where P, m (cos0) are normalized associated Legendre polynomi- 
als of degree / and order m, 0 is the colatitude, $ is the longi- 
tude, and r is the radius at which the field is evaluated. C/ m 
and Si,„ are the harmonic coefficients, GM =42828.28 km 3 s“ 2 
and R = 3394.2 km [Smith et al, 1993]. Because of the highly 
elliptical orbits and various inclinations of the Mariner 9 and 
Viking Orbiters, the actual resolution of the gravity field varies 
considerably with latitude and to a lesser extent with longitude. 
The resulting field has a formal half-wavelength horizontal 
resolution of 215 km, although in practice, the field only ap- 
proaches this resolution where low-altitude spacecraft tracking 
exists to constrain the field. The region of low-altitude cover- 
age extends roughly from 50 “N to 40 °S latitude, with the ef- 
fective resolution of the gravity field being considerably de- 
graded as one approaches the poles. Formal error estimates for 
the free-air gravity range from about 50 mGal in the equatorial 
zone to in excess of 80 mGal at the South Pole [Smith et al, 
1993]. Another spherical harmonic degree 50 model of the 
Mars gravity field, MarsSOc, has recently been developed 
[Konopliv and Sjogren, 1995]. GMM1 and Mars50c differ 
significantly from one another in their high-degree harmonics, 
but are quite similar in their low-degree harmonics. Because 
the emphasis in this paper is on the low-degree harmonics, we 
consider only GMM1; our inversion results are essentially un- 
changed if Mars50c is used instead. The free-air gravity ano- 
maly is 8g = - BU/dr . 


We also use a spherical harmonic degree 50 expansion 
of the topography of Mars [Bills and Nerem, 1995], of the form 


tf(G,<t>) = R 


50 / 

1+Z X P/,„(cos9) 

1 = I m = 0 


x (A/ m cosrnc(> + #/ m sin»i<t>) j (2) 

where H is the topography at a point and A )m and B im are the 
harmonic coefficients. The topography harmonics were derived 
from a digital version of the current Mars topography model 
[f/.S. Geological Survey (USGS), 1989]. The estimated vertical 
uncertainty in the input topography model varies from 1 km 
near the equator to more than 2 km near the poles [U5G5, 
1989]. The USGS topography model is referenced to a Ma- 
riner 9 era degree 4 equi potential surface [Jordan and Lorell, 
1975]. In contrast, the topography maps and harmonics in this 
paper are referenced to the degree 50 equipotential surface of 
Smith et at, [1993]. An additional difference between the 


USGS topography and the spherical harmonic model is the zero 
reference elevation. The USGS topography model is referenced 
to zero elevation at the 6.1-mbar pressure level, and in this 
reference frame the topography has a globally averaged eleva- 
tion of nearly 2 km. In contrast, in the spherical harmonic 
model the zero elevation is set by the mean radius R and any 
map based on the spherical harmonics will produce a surface 
with a globally averaged elevation of zero, regardless of the 
range of harmonic degrees used in the expansion. Thus numer- 
ical values from our topographic maps can be directly com- 
pared with the values on the USGS maps only if this offset is 
accounted for. We will sometimes use the shorthand notation 
U fm to refer collectively to a pair of coefficients C/„, and S fm 
and the notation H ]m to refer collectively to a pair of 
coefficients A tm and B tm at a particular value of / and m. 


Statistical Analysis of Harmonic Models 

In this section, we briefly discuss some of the statistical 
properties of the Martian gravity and topography harmonic 
models. We consider only the long-wavelength portions of 
these fields, emphasizing issues that are relevant to the inver- 
sion developed in the next section. Analysis of the full har- 
monic degree 50 fields was presented by Bills and Nerem 
[1995]. Comparisons of Mars with the gravity and topography 
of Earth, Venus, and the Moon have been discussed by Esposi- 
to et al [1992] and Balmino [1993]. 

The circles in Figure 1 show the gravitational power (the 
sum of the squares of the coefficients) as a function of harmon- 
ic degree. The triangles in Figure 1 show the topography 
power, expressed in terms of the gravitational potential due to 
uncompensated topography [Bills et al, 1987, equation (50)]. 
For degrees 3 to 10, the gravity power is typically a factor of 



Spherical Harmonic Degree 

Figure 1. Root-mean-square amplitude of the gravitational po- 
tential of Mars (circles) and of the gravitational potential due to 
uncompensated topography (triangles) as a function of spherical 
harmonic degree. 
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about 2 to 3 less than the uncompensated topography power, 
implying that the long-wavelength topography of Mars is com- 
pensated to a considerable degree. Between harmonic degrees 
14 and 19, the gravity power is always between 71% and 87% 
of the uncompensated topography power, implying that topog- 
raphy in this wave band is largely uncompensated. Maps of the 
geoid anomaly, free-air gravity anomaly, and topography in the 
degree 14 to 19 wave band are dominated by high-amplitude 
positive anomalies centered on Olympus Mons, Arsia Mons, 
Ascraeus Mons, Alba Patera, and Elysium Mons. These vol- 
canos are all partially supported by elastic flexure [Comer et 
ai, 1985 \ZuberetaL, 1993; Kiefer et ai, 1995]. 

A more quantitative way to analyze the relationship 
between gravity and topography is to consider a linear relation- 
ship of the form 


Utm = Ft Him + hm (3) 

Here, F t is the admittance and I /m are the coefficients for the 
part of the gravitational potential which is uncorrelated with the 
topography. The admittance may be estimated using standard 
least squares methods [e.g.. Bills et ai, 1987, appendix]. This 
is shown with one standard deviation uncertainties in Figure 2. 
For topography which is Airy compensated at a depth D below 
the surface, the admittance is given by [e.g.. Bills et ai, 1987] 


F f = 


3 Pc 


(2/+ 1 ) p a , 




(4) 


where p c is the assumed crustal density (2900 kg m" 3 ) and p av 
is the average planetary density (3930 kg m~ 3 ). This is shown 
in Figure 2 for D = 100 km (bottom dashed line), D = 200 km 
(middle dashed line), and £> = 400 km (top dashed line). Clear- 
ly, no single compensation depth can satisfy the observations 
[e.g., Phillips and Saunders, 1975; Lambeck, 1979]. The de- 
gree 2 and 3 admittances require compensation depths of 1400 
and 550 km, respectively, which vastly exceeds any reasonable 
estimate of the crust’s thickness and indicates that dynamic 
support by mantle convection may play an important role. 
Moreover, the admittance produced by a convecting system can 
decline much more rapidly as a function of harmonic degree 



Spherical Harmonic Degree 

Figure 2, Observed admittances as a function of spherical har- 
monic degree. One sigma uncertainties are also shown for de- 
grees 2 through 4. For higher degrees, the uncertainties are 
smaller than the plotted symbol size. The dashed lines are ad- 
mittance spectra for Airy compensation at depths of 100 km 
(bottom line), 200 km (middle), and 400 km (top line). 



Spherical Harmonic Degree 

Figure 3. Correlation between gravitational potential and to- 
pography on Mars as a function of spherical harmonic degree 
(circles and solid lines). Dashed line is 95% statistical 
confidence limit. 


than for either Airy or Pratt compensated layers [Kiefer et ai, 
1986; Kiefer and Hager ; 1991]. The rapid decline in admit- 
tance observed for degrees 2, 3, and 4 is therefore also qualita- 
tively consistent with an important role for dynamic support. 
From degrees 4 to 20, Airy compensation at a depth of 100 km 
can account for most of the observations. However, this does 
not rule out a role for convection at these wavelengths. As 
shown later in this paper, most of the gravity field and some of 
the topography may be supported by mantle convection, at least 
up to the /= 10 cut off used in our density inversions. 

The admittance values in Figure 2 are only meaningful if 
the linear relationship between gravity and topography (equa- 
tion (3)) is a good fit to the data. This can be assessed using 
the linear least squares correlation coefficient. In Figure 3, the 
circles and solid lines show the correlation between gravitation- 
al potential and topography as a function of harmonic degree, 
while the dashed line shows the 95% statistical confidence lim- 
it. The correlation is statistically significant for degrees 2, 5, 7, 
and 10 to 19 at the 95% confidence level and for degrees 3 and 
6 at the 90% confidence level. Although the correlation is not 
statistically significant at a few degrees, enough degrees are 
well-correlated that we can be sure that the general shape of the 
admittance spectrum (a rapid decline at low degrees and nearly 
constant admittance thereafter) is well constrained. The error 
bars on the admittances in Figure 2 also demonstrate this. 

Modeling Procedures 

In this section, we describe our method for inverting 
gravity and topography for models of the internal mass distribu- 
tion of Mars. Results from our inversions are presented in the 
following section. Given the gravitational potential of a planet, 
one can always solve for density anomalies at any specified 
depth in the planet which will exactly reproduce the observed 
potential. However, the resulting mass distribution will not in 
general be consistent with the planet’s topography and common 
concepts such as isostasy, flexure, or mantle convection. If 
density anomalies are allowed to exist in two shells at different 
depths, then it is possible to find a mass distribution which will 
satisfy both the observed gravitational potential and the topog- 
raphy. The results of such two layer models are non-unique. 
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depending on the assumed depths of the density anomalies and 
on the assumed compensation mechanisms. Nevertheless, by 
considering a range of plausible models, it is possible to deter- 
mine the range of density distributions which are consistent 
with the observed gravity and topography of a planet. 

Many prior quantitative analyses of the gravity and to- 
pography of Mars have focused on the role of crustal and 
lithospheric processes. Many of these studies have emphasized 
the Tharsis region [Sleep and Phillips, 1979, 1985; Banerdt el 
ai, 1982, 1992; Willemann and Turcotte, 1982; Finnerty et ai, 
1988; Phillips et ai, 1990]. Other studies have considered 
Elysium [Janie and Ropers, 1983], large shield volcanos 
[Sjogren, 1979; Zuber et ai, 1993; Kiefer el ai, 1995, Turtle 
and Melosh, 1995], the crustal dichotomy boundary [Janie, 
1983; Phillips, 1988], and impact basins and large craters 
[Sjogren and Wimberley, 1981; Sjogren and Ritke, 1982; Solo- 
mon et ai, 1983]. In several of the Tharsis studies, the possi- 
ble contribution of mantle convection to the gravity and topog- 
raphy was explicitly rejected, although none of the studies cited 
above actually attempted to model the effects of mantle convec- 
tion. However, the prolonged history of volcanistn in the 
Tharsis and Elysium regions [Plescia and Saunders , 1979; Scott 
and Tanaka, 1980; Tanaka and Scott, 1987; Plescia, 1990] can 
be readily explained by adiabatic decompression of upwelling 
material. This implies that mantle convection probably has 
been an important process, at least in these two areas of Mars, 
perhaps for much of the history of the planet. 

The density anomalies associated with this convection 
will contribute to the gravity and topography. If the elastic 
lithosphere is very thick, as is true for some of the models cited 
above, then there may be little dynamic topography produced 
by the convective flow. Nevertheless, the density anomalies 
within the convecting mantle will contribute to the gravitational 
potential; this effect has not been incorporated in any of the 
models cited above. At most, these models allowed for static, 
Pratt-compensated density anomalies at depths of up to a few 
hundred kilometers. However, at the lowest harmonic degrees, 
density anomalies below the middle of the mantle may still 
contribute significantly to the gravitational potential [ Richards 
and Hager, 1984]. In addition, the gravity-topography admit- 
tance for convecting systems is generally not the same as for a 
Pratt-compensated layer [Kiefer et ai, 1986; Kiefer and Hager, 
1991]. Existing models which neglect mantle convection there- 
fore should not be considered as complete or definitive models 
for the gravity and topography of Mars, but rather as explora- 
tions of one type of end-member model. 

In this same spirit of exploring end-member models, we 
consider here another plausible but previously unexplored type 
of end-member model, in which a combination of mantle con- 
vection and isostatically compensated crustal structure exactly 
satisfies the long-wavelength gravity and topography. For rea- 
sons explained in greater detail below, we do not include the 
effects of elastic flexure of the lithosphere in this study. In real- 
ity, it is likely that convection, isostasy, and flexure all play 
roles in producing the gravity and topography of Mars, and ul- 
timately a model which simultaneously incorporates all three 
processes will be necessary. However, such a model will of 
necessity be fairly complex, and it is useful to first understand 
how each of the end-member mechanisms may contribute 
separately to the overall picture. The models presented here 
describe one such end-member. In the longer term, acquisition 
of improved gravity and topography data sets [Smith et ai, 


1990; Zuber et ai, 1992] and of seismic and heal flow data 
from the surface may be necessary to fully resolve the relative 
roles of isostasy, flexure, and convection. In addition to the 
inversions presented here, several studies have developed finite 
element simulations of mantle convection on Mars [Kiefer and 
Hager, 1989; Schubert et ai, 1990]. The finite element models 
require fewer approximations but do not exactly reproduce the 
observed gravity and topography. Our inversions require a 
greater degree of parameterization but do produce density 
models which exactly reproduce both the observed gravity and 
topography. Moreover, the inversions are computationally 
much more efficient, allowing us to rapidly explore the effects 
of varying model parameters. 

Inversion Method 

Recently, several two-layer models for the gravity and 
topography of Venus have been presented in which the deeper 
density layer drives viscous mantle flow and the upper density 
layer represents crustal structure [Bills and Fischer, 1992; 

Grimm and Phillips, 1992; Herrick and Phillips, 1992]. The 
models described here also combine mantle flow and crustal 
structure. Our approach most closely resembles the global 

Venus model of Herrick and Phillips [1992], although our im- 
plementation differs in some details. We split the observed 
gravity and topography into two contributions: 

U? m +U[ m =U lm (5a) 

+ = (5b) 

Here, U Un and H {m are the observed gravitational potential and 
topography at degree / and order m; recall that this notation im- 
plicitly refers to both the sine and cosine terms at a given / and 
m. The terms with superscript d are contributions to the gravi- 
ty and topography from the deep density anomaly shell. The 
terms with superscript s are contributions to the gravity and to- 
pography from the shallow density anomaly shell. Because we 
seek to derive information about variations in mantle density 
and crustal thickness, equations (5)-(7) are all expressed in di- 
mensional terms. Thus the Smith et ai [1993] gravity harmon- 
ics are multiplied by GM/R prior to being inserted for U lm on 
the right side of equation (5a). Similarly, the Bills and Nerem 
[1995] topography harmonics are multiplied by R before being 
inserted for H /m on the right side of equation (5b). 

We assume that density anomalies in the deep layer 
represent mantle temperature variations associated with convec- 
tive flow. This flow deforms the surface, producing topography. 
The gravitational potential is the sum of contributions from the 
mantle mass anomalies and from mass anomalies associated 
with flow-induced topography at the surface and at the base of 
the mantle (all three contributions are summed together in the 
weighting function X { in equation (6a)). The gravitational po- 
tential and topography produced by the convective flow can be 
written as 

Ui!„ = J XAr) 5p f m (r) dr (6a) 

R , 

Hi' m = — jnr)dpf m (r)dr (6b) 

P.sur r 

In equations (6a) and (6b), G is the gravitational constant, R is 
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the planetary radius, R\ and R 2 are the radii of the top and bot- 
tom of the shell in which density anomalies are located, and 
p sur is the density of the uplifted material at the surface. Both 
the surface of density contrast p c and the crust-mantle interface 
of density contrast (pm-pj are uplifted. Both interfaces contri- 
bute to the uplifted mass anomaly, so the sum of the two densi- 
ty contrasts, the mantle density p m , is used as the value for p sur . 
The functions Xfr) and 7)(r) are weighting functions, which are 
derived by solving the continuity and conservation of momen- 
tum equations for the viscous flow problem. We assume a 
depth-dependent, Newtonian rheology and therefore solve for 
X } (r) and T ( (r) using a propagator matrix technique [Richards 
and Hager , 1984]. The propagator matrix formalism allows 
viscosity to vary with depth but not laterally. In the absence of 
lateral viscosity variations, X { and 7} are independent of spheri- 
cal harmonic order m . The values of X f and T } are functions of 
the assumed boundary conditions and of the depth-dependence 
of the viscosity. We assume a no-slip (zero horizontal veloci- 
ty) top boundary and a tree-slip (zero shear stress) boundary at 
the base of the mantle. The assumed mantle viscosity model is 
discussed in greater detail below. 

Recent two-layer inversions of Venus [Bills and Fischer , 
1992; Grimm and Phillips, 1992; Herrick and Phillips, 1992] 
have ail assumed loading at a specified depth in the mantle, in 
essence making 8p f m (r) a delta function. Herrick and Phillips 
did consider various choices of loading depth, but each of their 
model inversions used a single fixed loading depth. In contrast, 
we believe it is important to allow loading of the mantle to oc- 
cur over a range of depths in each inversion. This is because 
the various harmonic degrees have different depth sensitivities, 
and because for any given degree, the functions X/(r) and 7}(r) 
have different depth dependences [Richards and Hager, 1984]. 
We do not have sufficient constraints to independently deter- 
mine Sp at more than two depths, but we can specify an a 
priori functional form for 8p (r) and simply solve for a scaling 
constant. We follow Kiefer ei al. [1986] in assuming that the 
loading is independent of depth in the mantle, that is 8p(r) = 8p. 
This is a reasonable first approximation to some types of con- 
vective structures, such as mantle plumes. However, it is likely 
that there are vertical variations in the mantle thermal structure 
of Mars, so our results are best considered as a vertical aver- 
age. Although a more detailed model for 5p(r) would be desir- 
able, we think that the simplified form used here is preferable 
to ignoring mantle flow altogether, as has been done in all prior 
quantitative studies of Martian gravity and topography. 

Jarvis and Peltier [1986] showed that when the tempera- 
ture distributions at various depths in a convecting layer are ex- 
panded in the wavenumber domain, the thermal boundary 
layers at the top and bottom of the system have different spec- 
tra from the remainder of the system. This implies that the sim- 
ple approximation that 8p(;) is independent of depth can not be 
extended to the boundary layers. For this reason, we do not in- 
corporate the upper thermal boundary layer into the integration 
of equations (6a) and (6b). Instead, we treat the upper thermal 
boundary layer as part of our shallow density layer and set the 
integration limit R 2 to be the average base of the upper thermal 
boundary layer. The Earth’s thermal boundary layer reaches an 
asymptotic value of 95 km [Stein and Stein, 1992], with an 
average thickness of roughly two-thirds of this, or about 60 km. 
Scaling the heat flux by planetary mass and inversely by sur- 
face area suggests that Mars has a heat flux that is about 40% 
of Earth’s heat flux. The upper boundary layer thickness scales 


inversely with heat flux, so we adopt 150 km as the value for 
R 2 . Strictly speaking, we should also handle the lower thermal 
boundary layer density model separately from the overlying 
mantle. However, because both X/(r) and T f (r) go to zero at the 
base of the convecting layer, the lower thermal boundary layer 
makes a negligible contribution to the integrals in equations 
(6a) and (6b). Thus the details of how the lower thermal boun- 
dary layer is included in 8p(r) are unimportant. 

Our shallow density anomaly layer includes both density 
variations due to lateral variations in the structure of the upper 
thermal boundary layer and density variations due to variations 
in crustal thickness. Density anomalies in the thermal boundary 
layer create topography by viscous flow, and their gravitational 
and topographic effects can therefore be treated by the formal- 
ism in equations (6a) and (6b). We noted earlier that Pratt com- 
pensation is not a good analog for a complete convecting sys- 
tem. However, if one considers just the upper thermal boun- 
dary layer of a convecting system, where the loading depth is 
small in comparison with the horizontal wavelength, then the 
viscous flow calculation gives results which are essentially the 
same as for the case of Pratt compensation. In other words, 
within the upper thermal boundary layer, the function T f is 
close to its limiting value of -1,0. In this work, we also as- 
sume that variations in crustal thickness are compensated by 
Airy isostasy. For wavelengths that are much longer than the 
vertical thickness of the isostatic compensation layer, Airy and 
Pratt compensation are indistinguishable on the basis of their 
gravity signatures. For this reason, it is necessary to consider 
both boundary layer thermal structures and crustal thickness 
variations as part of the shallow layer. This view contrasts 
with recent Venus inversions cited earlier, where it was as- 
sumed that the effects of the crust and mantle could be separat- 
ed into different density layers. 

Although the foregoing discussion shows that both the 
upper thermal boundary layer and the crust may contribute to 
our shallow density layer, it is not possible to specify a priori 
the relative importance of these two sources. Indeed, the rela- 
tive strength of the two sources may be different for each term 
in the harmonic expansion. For the sake of formulating a 
mathematically well-posed problem, we therefore formulate our 
inversion assuming that the shallow layer is due solely to crus- 
tal thickness variations. However, when interpreting our results, 
we must remember that boundary layer temperature variations 
are also part of the shallow layer. When the shallow layer is 
formulated in terms of crustal thickness variations, it is con- 
venient to let the topography supported by the shallow layer be 
the independent variable. Thus we may write 


UL= 


47tGR 
21 + \ 


Pc ML 




(7a) 


HL = HL (7b) 

Equation (7a) is the relationship between gravitational potential 
and topography for a crust of density p c , assuming that the to- 
pography is Airy compensated at a depth D below the surface 
of the planet. The admittance defined in equation (4) is essen- 
tially a non-dimensional version of equation (7a). Equation 
(7b) simply emphasizes that topography is the independent vari- 
able for this layer of the inversion. By substituting equations 
(6a), (6b), (7a), and (7b) into equations (5a) and (5b), we end 
up with a linear system of two equations in two unknowns, 
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a) Mars Topography: degree 2-10 



b) Mars Gravity Anomalies: degree 2-10 



Figure 4. (a) Observed topography of Mars, filtered to include 
only harmonic degrees 2 to 10. The contour interval is 1 km, 
with lows shaded, (b) Observed free-air gravity anomalies of 
Mars for harmonic degrees 2 to 10. The contour interval is 100 
mGal, with lows shaded. The small circles in these and the 
other maps represent eight large shield volcanos: Arsia Mons, 
Ascraeus Mons, Olympus Mons, Pavonis Mons and Alba Patera 
in Tharsis and Albor Tholus, Elysium Mons, and Hecate Tholus 
in Elysium. 


which is easily inverted for 5p f m and H} m . This process is car- 
ried out separately for each term in the harmonic expansion. 

Because this is an exploratory, global-scale study, we 
have restricted our inversions to harmonic degrees 2 to 10, 
corresponding to a half- wavelength resolution of 1065 km. The 
observed gravity and topography of Mars filtered to this resolu- 
tion arc shown in Figure 4. Higher resolution studies of some 
features on Mars are certainly worth pursuing, but are probably 
best handled in the context of regional rather than global 
models. Given the resolution of our global model, we em- 
phasize the interpretation of large-scale geologic provinces: ma- 
jor impact basins, the crustal dichotomy, Tharsis, and Elysium. 

In this work, we do not consider the effects of elastic 
lithospheric flexure. This approximation is justified for several 
reasons. First, the effects of flexure are most important at short 
wavelengths, whereas we are considering the longest 
wavelengths on the planet. Second, the specific nature of some 
of the features of interest here suggests that flexure is unimpor- 
tant. For example, the crustal dichotomy clearly formed very 
early in the history of Mars. Theoretically, the mantle heat 
flow should have been high at that time, and the elastic litho- 
sphere should have been quite thin. Observationally, if flexure 


were important in supporting the crustal dichotomy, then the 3 
to 4 km of relief across this boundary would be accompanied 
by a free-air gravity anomaly of about 400 mGal. The absence 
of such a pronounced gravity anomaly along most of the dicho- 
tomy boundary indicates that flexure is unimportant in support- 
ing the dichotomy, at least at the wavelengths relevant to this 
study. The major impact basins also formed early in Martian 
history. Flexure should therefore be unimportant for the main 
basin topography, both because the heat flow was high and be- 
cause impact-induced fracturing further reduced the 
lithosphere’s strength. The long-wavelength gravity anomaly of 
only -60 mGals at the 7-km-deep Hellas basin is consistent 
with a lack of flexural support for basin topography. However, 
as the lithosphere cooled and thickened with time, it might pro- 
vide some flexural support for later loads emplaced within the 
basins. We discuss this below in the context of the Utopia 
basin mascon. 

Third, even if clastic flexure were at one time capable of 
supporting substantial vertical loads, viscous relaxation of elas- 
tic stresses will decrease the support of such loads as time 
progresses. The general trend of planetary thermal evolution is 
cooling of the interior and thickening of the elastic lithosphere 
with time. We noted earlier that the prolonged volcanic history 
at Tharsis and Elysium implies that convective upwelling has 
probably occurred in these regions for most of the last 4 billion 
years. In the first billion years of Martian history, the litho- 
sphere must have been quite thin [Schubert and Spohn . 1990] 
and incapable of resisting convective uplift of the surface. 
Later, as the planet gradually cooled and the Rayleigh number 
declined, the convective stresses and hence surface uplift gradu- 
ally increased [W.S. Kiefer, manuscript in preparation. 1996). 
The timescale for this cooling is about 10 9 years [Schubert and 
Spohn, 1990], Although the elastic lithosphere tends to resist 
this increasing convective uplift, the very slow timescale allows 
viscous relaxation to relieve most or all of the resisting elastic 
stresses in the lithosphere. Accordingly, surface uplift due to 
mantle convection is likely to remain in quasi-equilibrium with 
the underlying thermal structure. Because viscous relaxation 
occurs most rapidly at long wavelengths [e.g., Solomon et ai, 
1982], there is no inconsistency between this argument and the 
observation that elastic flexure contributes to the support of 
shorter-wavelcngth features such as shield volcanos [Comer et 
ai, 1985; Zuber et ai, 1993; Kiefer et ai, 1995]. The behavior 
of the gravity power spectrum (Figure I ) noted earlier is in fact 
quite consistent with little or no flexural support of topography 
at the lowest harmonic degrees and greater flexural support at 
high harmonics. 

Structural and Viscosity Models 

An important consideration in performing this inversion 
is the internal structure of Mars, particularly the thickness of 
the mantle. At present, the only geophysical constraints on the 
core and mantle structure of Mars come from the mean density 
and the moment of inertia. The mean density is well con- 
strained, but the moment of inertia is uncertain, with estimates 
ranging between dimensionless values of 0.345 [Bills, 1989] 
and 0.365 [Kcutla, 1979]. This uncertainty allows a broad 
range of possible internal structures. Even if the moment of in- 
ertia is specified exactly, uncertainty in the internal structure 
remains. This is because a two-layer (mantle and core) model 
requires specification of the core radius, the core density, and 
the mantle density, hut only two constraints (density and mo- 
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ment of inertia) are available. We consider a reference struc- 
tural model with a dimensionless moment of inertia of 0.355, 
midway between the values proposed by Kaula [1979] and Bills 
[1989]. This value of the moment of inertia was also recently 
advocated on geochemical grounds by Longhi el al. [1992]. 
Because the two-layer model is underconstrained, we choose to 
specify that the core’s radius is 1840 km, or 54% of the plane- 
tary radius, which is the same fraction of the planetary radius 
that the Earth’s core occupies. The core and mantle densities 
for this model are 7250 kg m -i and 3300 kg m“\ We consider 
only an average density for the entire mantle and neglect the 
small (“6%) density change across the olivine-spinel phase tran- 
sition. This has only a slight effect on the calculation of how 
the gravitational acceleration varies with depth and is unimpor- 
tant relative to the other uncertainties in the model. For this 
model, the hydrostatic contributions to the gravitational poten- 
tial are -7.585 x 10 -4 for C 2 « and 2.7 x KT 6 for C 4( , [Sleep and 
Phillips, 1985, Table B 2]. Hydrostatic contributions to the 
higher degree zonals are negligible. 

Another important consideration in this modeling is how 
the mantle’s viscosity varies with depth. We consider a refer- 
ence viscosity model that includes a high-viscosity near-surface 
layer (simulating the cold, high-viscosity thermal boundary 
layer) over an isoviscous mantle. The surface layer is 150 km 
thick and is 100 times the viscosity of the underlying mantle. 
The thickness of this layer is based on the heat flow scaling ar- 
guments outlined above. Finite element simulations have 
shown that once a high-viscosity lid is included in a convection 
model, an asymptotic state is reached such that increasing the 
thickness or viscosity contrast of the high-viscosity layer has a 
negligible effect on the geoid or topography produced by the 
convective flow [Kiefer and Hager, 1992]. Increasing the lid’s 
viscosity contrast from 100 to 10 4 in the present inversion has a 
negligible effect on the results. Increasing the lid thickness 
from 150 km to 250 km changes the mantle temperature range 
by less than 5% and the crustal thickness range by less than 
1%. These results indicate that our calculations arc not sensi- 
tive to the precise choice of high-viscosity lid parameters. The 


effects of alternative structural and mantle viscosity models are 
considered later. 

Results 

Reference Model 

Figure 5 shows the mantle temperature anomalies in- 
ferred from the inversion of the reference model. The inver- 
sion actually produces estimates of mantle density anomalies, 
which are converted to temperature anomalies using the thermal 
expansion relationship 

5p=-p m a87’. (8) 

Here, 5 p and 8 T are the density and temperature anomalies, p m 
is the unperturbed mantle density, and a is the thermal expan- 
sion coefficient (3 x I0“ 5o C H ). This relationship can be ap- 
plied separately to each harmonic coefficient. The results in 
Figure 5 are shown with elevated temperatures in white and 
cold regions shaded. The total range of temperature anomalies 
in Figure 5 is +350 K to -210 K, a reasonable range for a con- 
vening mantle. In comparison, an earlier treatment of Tharsis 
in terms of Pratt compensation would require lateral tempera- 
ture variations of 2000 K [Sleep and Phillips, 1979, 1985]. As 
Sleep and Phillips noted, such a large variation in temperature 
is implausible, with density differences due to compositional 
variations being more likely. In contrast with this earlier work, 
our model requires much smaller temperature variations be- 
cause these variations extend through the mantle rather than be- 
ing confined to the upper few hundred kilometers in a thermal 
boundary layer. 

The most prominent feature of Figure 5 is the pro- 
nounced region of elevated temperatures beneath Tharsis. 
There are two temperature maxima beneath Tharsis. One is 
just southeast of Pavonis Mons, but elongated in the directions 
of Arsia and Ascraeus Montes. The second is centered on 
Olympus Mons. The high-temperature region extends north- 
ward to include Alba Patera, but the temperature there is not as 


Mars Mantle Temperature: degree 2-10 



Figure 5. Temperature anomalies in the mantle of Mars for spherical harmonic degrees 2 to 10 for the reference 
model. The contour interval is 50 K, with negative anomalies shaded. 
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high as at either Olympus or Pavonis Montes. The existence of 
high mantle temperatures in Tharsis is quite consistent with the 
prolonged history of volcanism in this region. Indeed, the 
Tharsis region has been volcanically active in geologically re- 
cent times [Plescia and Saunders, 1979; Scott and Tanaka, 
1980] and may still be active. There is also a hot anomaly in 
the Elysium region, centered very close to Elysium Mons. The 
amplitude of this anomaly is smaller than the Tharsis anomaly, 
consistent with the lower level of volcanic activity in Elysium. 
The pattern of cold anomalies in Figure 5 does not generally 
correlate with the surface geology. 

The most surprising feature of Figure 5 is the prominent 
hot anomaly in the Utopia region. This correlates closely with 
the Utopia impact basin mapped by McGill [1989], The hot 
anomaly cannot be a relic of the impact formation process be- 
cause thermal diffusion will cool impact heating to negligible 
levels in just a few hundred million years [Bratt et ai, 1985], 
and the Utopia basin is likely to be nearly 4 billion years old. 
There is evidence for volcanic activity in the early Amazonian 
in this region, but it is believed due to flows emanating from 
Elysium Mons rather than being indigenous to the Utopia re- 
gion (unit Ael3 of Greeley and Guest [1987]). It is possible, 
however, that earlier Utopia-centered volcanism was covered 
over by the later Elysium eruptions. McGill [1989] noted the 
existence of a gravity anomaly which closely correlates with 
the position of the Utopia basin and suggested that this anoma- 
ly is a mascon. 

The hot anomaly in Utopia in Figure 5 may be the sig- 
nature of this mascon. In order to produce a free-air gravity 
high in a topographic low such as an impact basin, there must 
be material present which is at least partially uncompensated. 
Our model requires that shallow loads be isostatically compen- 
sated, which means that shallow loads cannot reproduce the ob- 
served gravity. In order to produce the observed gravity high, 
the inversion therefore places a load in the mantle. With the 
viscosity model used here, a low-density anomaly is needed to 
produce a positive gravity anomaly [Richards and Hager , 
1984]. Thus a flexurally supported, shallow mass excess could 
show up in our inversion as a high-temperature anomaly in the 
mantle. In principle, one might also apply this argument to the 
high-temperature regions in Tharsis, and as noted earlier, 
several investigators have proposed flexural models for Tharsis. 
An important difference is that high mantle temperatures in 
Tharsis are consistent with (and possibly demanded by) the re- 
cent volcanic history, an argument that cannot be applied to 
Utopia. One possible difficulty with the mascon interpretation 
for Utopia is that proposed mascons in Hellas and Isidis 
[Sjogren and Wimberley , 1981; Solomon et ai, 1983] do not 
show the same signature in Figure 5 as the Utopia feature. 
This is not a resolution effect because both Hellas and Isidis 
show up clearly in the topography map at the same resolution 
(Figure 4a). However, the gravity anomalies in Hellas and 
Isidis are smaller in amplitude than in Utopia, so this may not 
be a serious objection. 

The observed free-air gravity anomaly in Utopia is about 
150 to 200 mGals relative to its surroundings [Smith et ai, 
1993]. By measuring the gravity anomaly relative to the sur- 
rounding region, we are implicitly assuming that there was no 
gravity anomaly prior to emplacement of the mascon load. 
This is a reasonable first-order approximation if the original im- 
pact basin was isostatically compensated, as argued earlier. 


The minimum load required to produce the mascon can be 
derived by assuming that the load is emplaced at the surface 
(either by volcanism or by sedimentary deposition) and is en- 
tirely uncompensated. In this case, the free-air gravity anomaly 
associated with the load is simply 

8 gU) = 2 Jt G pi, >a j 8 h(x), (9) 

where p )iud and 8 h(x) are the density and thickness of the load. 
This assumes a planar geometry for the load, but spherical 
geometry will not significantly alter this for a gravity anomaly 
which is 2000 km across. Assuming p lujd = 2900kg m -3 
(basalt), the load would need to be 1.2 to 1,6 km thick to pro- 
duce the observed gravity anomaly. The Elysium flows 
mapped by Greeley and Guest [1987] in this region could be a 
major contributor to this load. Sedimentary deposits surround 
the volcanic flows [Lucchitta et ai, 1986] and might also un- 
derlie the Elysium volcanic flows, thus contributing to the in- 
ferred load thickness. If the load density is less than assumed 
(e.g., unconsolidated sediments) or the load is partially compen- 
sated, then a greater load thickness is necessary to reproduce 
the observed gravity anomaly. Based on geological relation- 
ships in Utopia, Frey and Schultz [1990] estimated a minimum 
volcanic fill of 1 km, with a more likely value of 2 km. This is 
in good agreement with the estimate derived here, particularly 
if the load is partially compensated, 

Figure 6a shows the shallow layer loads for this inver- 
sion. The loads are shown in terms of isostatically supported 
surface topography, with a range of -7.7 to +4.8 km. The total 
variation in crustal thickness (surface relief plus relief on the 
crust-mantle interface) will be larger by a factor p m /(p m “p c ), 
where p c and p m are the crust and mantle densities. Pratt com- 
pensated thermal loads in the upper thermal boundary layer can 
also contribute to this shell. One can convert between Airy and 
Pratt supported loads using the relationship 

p c 5/r v = Lp m otSr v . (10) 

Here, 8/r v is the topography supported by Airy isostasy and 8T V 
is the boundary layer temperature anomaly averaged through a 
thermal lithosphere of thickness L. The superscript .v em- 
phasizes that these are the loads from the shallow layer portion 
of the inversion. Equation (10) allows one to convert from 
pure Airy compensation (Figure 6a) to pure Pratt compensation; 
in the more general case, a combination of both Airy and Pratt 
isostasy may apply at any given point. Crustal thickening is 
shown in white in Figure 6a and crustal thinning is shaded. As 
expected, the major impact basins (Hellas, Argyre, Isidis, Uto- 
pia, and Chryse) are all regions of thinned crust in Figure 6a. 
Except where affected by these basins, areas on the south side 
of the crustal dichotomy boundary are typically areas of thick- 
ened crust. Regions near the north pole also show up as thick 
crust, but this simply reflects a degree 1 component in the to- 
pography that is not included in our inversion. The region 
around Valles Marineris is thickened crust (note that the actual 
trough system cannot be resolved at spherical harmonic degree 
10). This crustal thickening is elongated to the east and west, 
in the direction of Valles Marineris itself. 

Figure 6a also shows considerable crustal thinning in 
Tharsis. The minimum shallow layer topography in the Tharsis 
region in Figure 6a is -5 km at Olympus Mons, implying a to- 
tal crustal thinning of about 41 km for crust and mantle densi- 
ties of 2900 and 3300 kg nT 3 . The shallow layer topography in 
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a) Mars Shallow Layer Topography: degree 2-10 



b) Mars Deep Layer Topography: degree 2-10 



Figure 6. Surface topography produced by the two density anomaly shells for spherical harmonic degrees 2 to 
10 for the reference model. The contour interval is 1 km, with negative topography shaded, (a) Topography 
supported by shallow shell, (b) Topography supported by deep shell. 


central Tharsis is only about -2 km, implying about 16 km of 
crustal thinning. Given the hot mantle temperature beneath 
Tharsis (Figure 5), we would also expect the thermal boundary 
layer to be hotter than normal. A hot boundary layer affects 
the gravity and topography in a manner analogous to thickened 
crust (equation (10)). Thus, if a hot boundary layer is present 
in Tharsis, the actual amount of crustal thinning must be even 
larger than shown in Figure 6a in order to compensate for the 
effects of the boundary layer. 

Considering the prolonged history of volcanism in 
Tharsis, the predicted crustal thinning seems unexpected. The 
most likely explanation is that mantle shear stresses at the base 
of the crust are able to thin the crust more rapidly than volcanic 
processes are able to thicken the crust. The prediction of crus- 
tal thinning in Tharsis resembles the model of Sleep and Phil- 


lips [1979]. Since the work of Sleep and Phillips [1979], crus- 
tal thickening in Tharsis has been generally preferred [Solomon 
and Head , 1982; F inner ty et ai, 1988; Phillips et ai , 1990]. 
These models all presume that the crust and low-density residu- 
um produced by magmatic activity remain confined in Tharsis. 
Our results suggest that this presumption may be incorrect. 
However, more detailed modeling is clearly necessary to assess 
the conditions under which basal shearing can thin the crust 
more rapidly than magma production can thicken it. Such 
modeling will require the addition of a magma production 
model to a finite element simulation of the temperature and 
velocity fields in the mantle. 

There are two ways in which the predicted crustal thin- 
ning in Tharsis could be reduced. First, if the regional topogra- 
phy in Tharsis is larger than assumed here (a possibility sug- 
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gested by recent interpretation of Viking radio occultation ob- 
servations [Smith and Zuber, 1995]), the actual and model 
dynamic topography would agree more closely, reducing the 
need for crustal thinning. Second, the mantle loading model 
assumes that the density anomalies are uniformly distributed 
with depth. Because the low-degree geoid weighting functions 
( X ( {r )) peak in the middle portion of the mantle and the 
corresponding topography weighting functions (T/(r)) always 
peak at the surface [Richards and Hager, 1984], changing the 
distribution of the density loading in equations (6a) and (6b) 
will alter our results to some extent. In particular, if there is 
enhanced loading in the middle of the mantle, the same geoid 
anomaly could be achieved with a lesser amount of dynamic to- 
pographic uplift, which in turn would reduce the required 
amount of crustal thinning in Tharsis. Such a loading model 
might arise in time-dependent convection, for example, from an 
instability which forms in the lower thermal boundary layer and 
rises toward the surface. Clearly, a range of models can satisfy 
the observed geoid and topography. Observations from a se- 
ismic network [Solomon et ai, 1991] ultimately may be re- 
quired to definitively determine the extent to which crustal 
thickening or thinning has occurred in Tharsis. 

Figure 6a shows that the shallow-layer topography has a 
small amplitude at the Elysium volcanos, about I km at Elysi- 
um Mons. The relatively small shallow-layer topography 
shown for Elysium could be entirely due to elevated thermal 
boundary layer temperatures, allowing the possibility that Elysi- 
um, like Tharsis, is a region of crustal thinning. However, the 
amplitude of crustal thinning must be less in Elysium than in 
Tharsis. Based on the observed topography and gravity (Figure 
4), the horizontal scale for convective flow beneath Elysium is 
smaller than for convective flow beneath Tharsis. Such differ- 
ences in the horizontal scale of convection have a substantial 
effect on the geoid and dynamic topography produced by the 
flow [Kiefer and Hager , 1992], and we speculate that it also 
may affect the relative balance of crustal thickening or thinning 
produced by the flow. The maximum shallow-layer topography 
amplitude in this area is located about 1000 km to the northeast 
of Elysium Mons. The amplitude of this maximum (3 km) is 
sufficiently large that for plausible values of the boundary layer 
temperature contrast (57"), some crustal thickening is required. 

Figure 6b shows the surface topography supported by the 
mantle layer. The planform is quite similar to the temperature 


Mars Shallow Layer Gravity Anomalies: degree 2-10 



Figure 7. Free-air gravity anomalies produced by the shallow 
shell for spherical harmonic degrees 2 to 10 of the reference 
model. The contour interval is 25 mGal, and negative 
anomalies are shaded. 
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Figure 8. Correlation coefficient as a function of harmonic de- 
gree between the topography supported by the shallow density 
layer and the topography supported by the deep density layer 
for the reference model. Dashed line is 95% confidence limit. 


anomalies in Figure 5, with pronounced topographic uplift in 
Tharsis. Comparing Figures 6a and 6b, it is clear that the two 
density shells support comparable amounts of total topographic 
relief. However, as discussed below, the relative contributions 
of the two shells to the topography is a strong function of 
wavelength. In contrast with the topography, the mantle shell 
is the dominant contributor to the observed gravity. The gravi- 
ty anomaly due to the deep shell is not shown because it is vir- 
tually identical to the observed gravity (Figure 4b). Figure 7 
shows the gravity anomaly due to the shallow shell, which has 
an amplitude of less than 25 mGal in most places, with the not- 
able exception of Hellas, where the shallow-shell anomaly (-80 
mGal) represents most of the total anomaly. The dominance of 
the deep shell in the gravity field is not surprising because the 
crustal loads are isostatically compensated at shallow depths, 
implying small associated gravity anomalies. 

Figure 8 shows the correlation between the topography 
supported by the deep and shallow shells as a function of har- 
monic degree. The correlation for the gravity produced by the 
deep and shallow shells is identical. The dashed line is the 95% 
statistical confidence limit for each harmonic degree. The 
correlation between the two shells is negative at all harmonics 
except degree 5, but only degrees 8, 9, and 10 are statistically 
significant at the 95% confidence level. In addition, degrees 2 
and 4 are significant at the 90% level. Recall that the shallow- 
layer shell may include contributions both from the upper ther- 
mal boundary layer and from crustal thickness variations. The 
thermal structure of the upper boundary layer should be posi- 
tively correlated with the deep-layer thermal structure. The gen- 
erally negative correlations shown in Figure 8 therefore are evi- 
dence for a significant crustal contribution to the shallow-layer 
density structure. 

Figure 9 considers the relative contribution of the two 
shells to the gravitational potential and the topography as a 
function of harmonic degree. This is expressed in terms of the 
fractional power in the gravity field due to the deep shell, 

X (ufj 2 

f(i )= — — — . (ii) 

ZOJfJ + WfJ 2 

m=() 
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Figure 9. Fractional contribution of deep density layer to the 
gravitational potential (circles) and the topography (triangles) of 
the reference model. 


Here, the numerator is the power in the deep-shell geoid and 
the denominator is the power of both the deep- and shallow- 
layer geoids. The fractional topographic power due to the deep 
shell can be defined in a similar manner. The circles in Figure 
9 are the deep layer fractional power for the geoid, which 
shows that the deep layer dominates the gravitational potential, 
with fractional power exceeding 90% at all wavelengths. On 
the other hand, the topography results (triangles) show that the 
lower shell strongly dominates only the low-degree topography, 
with the shallow shell dominating at /> 4. 

Variation of Model Parameters 

We have examined the effects of varying the planetary 
structural model, the mantle viscosity model, and the depth of 
the shallow-layer density anomalies on our model results. The 
map patterns produced by these various models are generally 
quite similar to those shown in Figures 5 to 7, with the results 
differing primarily in the amplitude of the required density 
anomalies. For the sake of brevity, we therefore omit the 
presentation of map results for these models and instead just 
summarize the amplitude results. 

Schubert and Spohn [1990] tabulated structural models 
of Mars for various choices of moment of inertia and core den- 
sity. To test the sensitivity of our choice of reference model 
parameters, we have considered their models with the largest 
and smallest core radii with our reference viscosity model. 
Their model with the smallest core has a core radius of 1546 
km, a mantle density of 3500 kg nf\ and a moment of inertia 
of 0.365. The range of mantle temperature anomalies in this 
model is -150 to +280 K, somewhat less than for our reference 
model. Decreasing the core radius below our reference value 
has a relatively small effect on the required density anomalies. 
This is because density anomalies below 1500 km depth in the 
Mars mantle have little affect on either the geoid or the surface 
topography. Thus the required range of &T in the mantle does 
not differ greatly in these two models. The difference is due in 
part to differences in the mantle density in the two models. The 
gravity and topography which the models must support are sen- 
sitive to the density variations in the mantle, but we have con- 
verted this to temperature anomalies here. This conversion is 
inversely proportional to the nominal mantle density (equation 
(8)), which therefore requires a larger range of 6 T in the refer- 


ence model, which has a smaller mantle density. Much of the 
remaining difference is due to differences in the hydrostatic 
correction to the C 2( » gravity term for the two different structur- 
al models. The range of crustal topography in the two models 
differs by less than 10%. Recently, Fei et al. [1995] have re- 
ported new mineral physics results that allow the core of Mars 
to have a radius of just 1300 km. Based on the similarity of 
the two model calculations described here, we do not think that 
further decreasing the core size will significantly alter our 
model results. 

The model with the largest core in the Schubert a/ul 
Spohn [1990] tabulation has a core radius of 2468 km, a mantle 
density of 2785 kg nrr\ and a moment of inertia of 0.345. 
This large core model requires mantle temperature anomalies 
ranging between -570 and +850 K, much larger than for the 
reference model. Because the mantle density in the large core 
model is 15% less than in the reference model, the conversion 
from density anomaly to temperature anomaly (equation (8)) 
implies larger temperature anomalies in the large core model. 
Also, the mantle is only 900 km thick in the large core model 
and 1550 km thick in the reference model. A larger tempera- 
ture contrast is needed to compensate for the reduced depth 
range of the integration (equations (6a) and (6b)) in the large 
core model. The Earth’s upper mantle temperature is about 
1700 K [Stein and Stein, 1992], and because of the strong 
temperature-dependence of rheology, if convective tlow is 
significant in the mantle of Mars, then the mean mantle tem- 
perature on Mars cannot be more than 100 to 200 K less than 
the Earth’s temperature. With temperature anomalies of up to 
+850 K relative to the mean mantle temperature, the large core 
model would require large-scale melting in the mantle, for 
which there is no evidence in the recent volcanic history of 
Mars. Accordingly, we must conclude either that such a large 
core does not exist on Mars or that mantle convection is not a 
primary contributor to the present-day long-wavelength geoid 
and topography of Mars. 

We have also considered the effect of varying the mantle 
viscosity structure. The isoviscous mantle used in the reference 
model is Venus-like [Kiefer and Hager, 1991], On the other 
hand, models of the Earth’s mantle viscosity structure increase 
strongly with depth [e.g., Hager and Richards, 1989]. As one 
test of the effects of a viscosity increase with depth on our 
results, we have used our reference structural model with a 
lower mantle that is a factor of 30 more viscous than the upper 
mantle, with the transition assumed to occur at a depth of 800 
km, near the middle of our reference mantle. Increasing the 
viscosity with depth in the mantle weakens the coupling 
between convective stresses and the surface. This decreases the 
amplitude of the dynamic topography produced by the convec- 
tive flow and also has a large effect on the geoid [Richards and 
Hager, 1984]. Including a low- viscosity upper mantle layer re- 
quires mantle temperature anomalies between -700 and +1300 
K. The large positive temperature anomalies required in this 
model are implausible, as was the case for the model with a 
large core. This implies either that Mars does not have a strong 
increase in mantle viscosity with depth or that mantle convec- 
tion is not the dominant contributor to the present-day geoid 
and topography on Mars. 

The compensation depth in equation (7a) represents a 
weighted average of the mean crustal thickness and the middle 
of the upper thermal boundary layer. In the reference model, 
this was assumed to be 50 km, but we have considered values 
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of this parameter between 25 and 100 km. Compensating the 
shallow-density anomalies at 25 km represents a situation 
where the shallow layer is predominantly due to crustal thick- 
ness variations in a crust with a small mean thickness. Com- 
pensating the density anomalies at 100 km could represent ei- 
ther a very thick crust or a situation in which the shallow-layer 
density anomalies primarily represent boundary layer tempera- 
ture variations. The inversion results are relatively insensitive 
to the particular value of D. Mantle temperature anomalies 
change by only 10% and the crustal topography changes by 
20% over this range of D. Increasing the compensation depth 
of the shallow shell increases the contribution of the shallow 
shell to the gravity field. However, even for £>=100 km, the 
deep shell dominates the gravity, with the shallow shell contri- 
buting no more than 30% of the power. The fractional contri- 
bution of the two layers to the topography is virtually indepen- 
dent of the assumed value of D. 

Tectonic Implications 

Another method for testing internal structure models of 
the sort described in this paper is to compare the surface 
stresses predicted by the models with the observed orientation 
of tectonic features. Because tectonic features in old terrains 
on Mars may be unrelated to the present-day gravity and topog- 
raphy, such comparisons have been confined to the Tharsis and 
Elysium regions [Banerdt et ai, 1982, 1992; Willemann and 
Turcotte, 1982; Sleep and Phillips, 1985; Hall el ai, 1986]. 
Some of the predicted stress patterns in these studies may be 
modified by considerations of the effects of time-dependent 
lithospheric loading [McGovern and Solomon , 1993] and of im- 
proved lithospheric failure criteria [Schultz and Zuher, 1994]. 
Although we have not made a detailed analysis of the horizon- 
tal stresses implied by the models in this paper, the following 
analysis suggests that our models are generally consistent with 
the state of stress implied by the observed tectonics. 

Because Tharsis is approximately axisymmetric [Wil- 
lemann and Turcotte, 1982], we have examined the surface 
stresses predicted by several spherical axisymmetric mantle 
convection simulations [ Kiefer , 1993; W.S. Kiefer and L.H. 
Kellogg, manuscript in preparation, 1996]. Although the details 
are complicated due to the time-dependent nature of these con- 
vection simulations, in general, the predicted stresses imply ex- 
tensional features (graben) oriented radial to Tharsis near the 
center of the uplift and compressional features (ridges) that are 
oriented concentric to Tharsis at larger distances from the 
center of uplift. This is generally consistent with the pattern of 
faulting observed in Tharsis. Inclusion of the stresses associat- 
ed with the predicted crustal thinning in Tharsis in our model 
will modify the predicted tectonics. However, because the to- 
pography produced by our mantle load is up to 4 times as large 
as the shallow-layer topography in central Tharsis, we expect 
that the crustal layer will only be a small perturbation on the 
tectonics predicted using just the mantle flow model. A more 
detailed treatment, including a discussion of how the time- 
dependent convective stresses may be related to the temporal 
history of faulting in Tharsis, will be presented elsewhere [W.S. 
Kiefer, manuscript in preparation, 1996). 

Conclusions 

We have presented a series of inversions of the gravity 
and topography of Mars for mantle and crustal structure. Our 


models incorporate the effects of viscous mantle flow (which 
was neglected in prior models of the long-wavelength gravity 
field ot Mars) and of isostatically compensated density 
anomalies in the crust and lithosphere. These models predict 
hot upwellings beneath Tharsis and Elysium, consistent with 
the long volcanic histories of these areas. Crustal thinning is 
predicted beneath Tharsis, which if correct, implies that thin- 
ning of the crust by mantle shear tractions has dominated over 
crustal thickening due to volcanism. The large impact basins 
are all regions ot thinned crust, and Utopia contains a mascon 
whose gravity anomaly is consistent with 1.2 to 1.6 km of un- 
compensated basalt, consistent with independent geological esti- 
mates of the amount of volcanic fill in this region. 

We have tested our model for sensitivity to assumed 
crustal thickness, mantle viscosity model, and core radius. The 
general planform of mantle and crustal structure does not 
depend on the assumed values of these parameters, but the am- 
plitudes of the required temperature and crustal thickness varia- 
tions do depend on these parameters. Models with a small to 
intermediate size core and without a low-viscosity zone in the 
upper mantle require lateral temperature variations of a few 
hundred degrees Kelvin. Models with a very large core or with 
an upper mantle low-viscosity layer require unreasonably large 
lateral variations in mantle temperature and can therefore be re- 
jected. These results are fairly insensitive to the mean crustal 
thicknesses and to the thickness and viscosity contrast of the 
near-surface layer, in all of these models, the mantle thermal 
structure is the dominant contributor to the gravitational poten- 
tial. The mantle also dominates the topography at the longest 
wavelengths, but harmonic degrees >4 are dominated by the 
crustal structure. The predicted stress distribution is generally 
consistent with the observed distribution of extensional and 
compressional features in the Tharsis region. 
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